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ABSTRACT 


This research article presents a comprehensive analysis of the distribution 
patterns of Rhododendron arboreum in Nepal, encompassing the current 
distribution and future projections. We employed an ensemble approach by 
combining eight selected models with AUC > 0.8. Resulting ensemble model 
demonstrated remarkable predictive accuracy, achieving an AUC value of 0.97 
and a TSS value of 0.77. Our analysis unveiled that within total land area of 
Nepal, approximately 115,317.1 km2 were deemed unsuitable for Rhododendron 
arboreum, while 32,198.89 km2 were identified as suitable habitats. Future 
projections for 2050 and 2070 under the SSP2-4.5 and SSP5-8.5 scenarios showed 
significant changes in habitat suitability. By 2050, under SSP2-4.5, 10,710.39 km2 
were projected to be lost, while 4,792.78 km2 were gained. Under SSP5-8.5, more 
significant loss of 14,237.64 km2 and a smaller gain of 5,696.41 km2 were 
anticipated. In 2070, further losses in suitable habitat were projected under both 
scenarios. A consistent upward movement in the suitable habitat elevation range 
was observed from the current period to the projected years of 2050 and 2070. 
These findings underscore the urgency of implementing proactive conservation 
measures to safeguard Rhododendron arboreum in the face of climate change. The 
research provides valuable insights for informing conservation strategies and 
emphasizes the ecological significance of this species in Nepal. Considering 
distribution patterns and elevation shifts, our study contributes to a 
comprehensive understanding of the species dynamics and supports effective 


conservation planning. 
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1. INTRODUCTION 


Over the last century, global temperatures have been discernible increased by 
approximately 1°C, with a particularly accelerated rise observed in the past three 
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decades (Zhao et al., 2021). This significant shift in climate patterns has profoundly impacted the spatial distribution of plant 
species. The dynamic nature of climate change has influenced and continues to shape the geographical distribution of plants on a 
global scale. The current magnitude of climate change poses such a significant threat that even the most stringent measures 
implemented to mitigate its detrimental effects would prove insufficient (Chhetri et al., 2021; Pachauri et al., 2014). The global 
climate shift poses a significant threat to Earth's biological diversity. As a result, the global climate shift has the potential to 
significantly impact the species distribution, leading to the loss of suitable habitats for some of them. Given the potential threats to 
species distribution and habitats, it is crucial to make informed management decisions. 

This requires a deeper understanding of how global climate change affects species distribution. These decisions are essential to 
minimize the destructive consequences of future climate change on species habitats (Ardestani and Ghahfarrokhi, 2021; Pressey et 
al., 2007). Species distribution models (SDMs) have become a robust instrument for investigating the geographic distribution of 
diverse life forms, including plants and animals. They have proven helpful in assessing the conservation implications of 
distribution patterns and predicting future climate change impacts. SDMs offer novel perspectives for evaluating species range 
shifts under different climatic conditions between present and future periods (Chhetri et al., 2021; Ardestani and Ghahfarrokhi, 
2021). These models utilize quantitative equations that incorporate environmental variables, species presence data, and, in some 
cases, absence to forecast species distribution (Crego et al., 2022). 

Rhododendron arboreum, commonly known as the Rhododendron or Laligurans in Nepal, is a species of flowering plant that holds 
immense cultural and ecological significance in the country. It is Nepal's national flower and deeply embedded in local traditions 
and festivities. It is one of the most stately and impressive Rhododendron species (Srivastava, 2012). The species can be observed in 
various regions of Nepal, including the western, central, and eastern areas, specifically in hilly open areas at altitudes ranging from 
1300 to 3200 meters. Its distribution also extends beyond Nepal to other parts of the Himalayas, spanning from Kashmir to 
Arunachal Pradesh (Shrestha, 2016; Srivastava, 2012). It is an evergreen tree that exhibits a dense branching pattern, reaching a 
height of up to 14 meters with a circumference of 2.4 meters. It flowers from March to April or June to September, displaying 
vibrant flowers ranging from deep red or crimson to pale pink. 

Reports indicate that the bright red forms of Rhododendron arboreum are commonly observed at lower elevations (Orwa et al., 
2009). This plant thrives in light (sandy) to medium (loamy) soil and requires a moderately moist and acidic soil environment 
(Srivastava, 2012). The wood grain of Rhododendron arboreum is utilized in the production of handles for the traditional Nepalese 
knife called ‘khukuri', packsaddles, gift boxes, gunstocks, and posts. Its durable and versatile wood is prized for its suitability in 
crafting these items (Paul et al., 2005). This species stands out as a favored and extensively harvested species within the 
rhododendron genus, primarily due to its exceptional qualities as fuelwood as it boasts a remarkable combination of high calorific 
value and an elevated fuelwood index value (Ranjitkar et al., 2014). In hilly regions, the flowers of Rhododendron arboreum, known 
for their sweet and sour taste, are utilized in making squash, jams, jellies, and local brew. Its flowers are a popular and delightful 
beverage, consumed daily as a refreshing appetizer and for their potential to prevent high altitude sickness (Pant et al., 2009; 
Prakash et al., 2007; Singh and Chatterjee, 2022; Solanki et al., 2013). 

Despite the well-established impact of climate change on various plant species, understanding its effects on Rhododendron 
arboreum in Nepal remains limited. This knowledge gap can be attributed to several factors. Previous research on climate change 
impacts has often focused on more widespread or economically significant plant species, giving relatively less attention to 
Rhododendron arboreum. Additionally, the species distribution range, which includes diverse geographical and altitudinal gradients 
within Nepal, presents challenges for comprehensive studies. Moreover, the cultural and ecological significance of Rhododendron 
arboreum in Nepal adds a unique dimension to its analysis, warranting specific research to assess the potential implications of 
climate change in its distribution and conservation. 

Given these reasons, there is a pressing need to fill this knowledge gap and generate insights into the effects of climate change 
on Rhododendron arboreum in Nepal. Predicting species distribution changes and their responses to climate change can offer valuable 
insights for ecological monitoring, species conservation, and effective management strategies (Bahadur et al., 2023; Moradi et al., 
2019). The main objective of this research article is to investigate the distribution patterns of Rhododendron arboreum in Nepal under 
current and future climate change scenarios. By assessing the potential impact of climate change on the distribution of this species, 
our study aims to provide valuable insights into the possible range shifts and conservation implications for Rhododendron arboreum 
in Nepal's changing climate 


Species 24, e81s1590 (2023) 2 of 14 


ANALYSIS ARTICLE OPEN ACCESS 


2. MATERIALS AND METHODS 

Study Area 

Situated in South Asia, Nepal is a mountainous nation spanning an area of 147,516 km2 and positioned between the latitudes and 
longitudes of 26°22’-30°27' N and 80°04’-88°12' E respectively (Figure 1) (Dhami et al., 2023). Nepal boasts a rich biodiversity that 
flourishes across its steep altitudinal gradient, ranging from 60 to 8,848 meters above sea level, owing to its diverse climate and 
topography (Bhattacharjee et al., 2017). This variation in elevation creates distinct habitats across five primary physiographic 
regions: the High Himalaya, High Mountain, Middle Mountain, Siwalik, and Terai (Uddin et al., 2015). The area experiences a mild 
climate characterized by dry winters and rainy summers (Karki et al., 2016). The region receives an average annual precipitation of 
1768 mm, accompanied by mean annual temperatures of 18 °C (Shrestha et al., 2000). 
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Figure 1 Study area (Nepal) with occurrence points of Rhododendron arboretum 


Data Collection 

Species Occurrence Data 

The acquisition of species occurrence data (n=232) relied on detecting and surveying established permanent sample plots during 
regular inventories conducted by the Forest Research and Training Centre (FRTC) under the Ministry of Forest and Environment in 
Nepal. This approach ensured the reliability and representativeness of the dataset for subsequent analyses and research on species 
distributions, thanks to the use of standardized protocols, taxonomic expertise, and repeated surveys over time. Additionally, we 
augmented the dataset with 71 occurrence records obtained from the Global Biodiversity Information Facilities (GBIF, 2023). 

Thus, a total of 303 occurrences of Rhododendron arboreum were collected. We manually eliminated occurrences that coincided 
with buildings and water bodies to refine the dataset. Additionally, to address spatial autocorrelation, we employed the ‘fuzzySim’ 
package in R Barbosa, (2015), to randomly select one occurrence per pixel of the predictor variables. As a result, our dataset 
consisted of 160 occurrence point that exclusively represented the presence of the species. 
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Current Environmental Parameter 

For the assessment of the current distribution of Rhododendron arboreum, a total of 23 environmental parameters were selected. These 
parameters comprised 19 bioclimatic variables and three topographic variables (elevation, slope, and aspect). The bioclimatic 
variables and the Shuttle Radar Topographic Mission (SRTM) Digital Elevation Model (DEM) were retrieved from WorldClim 
(http://www.worldclim.org/) at a resolution of 30 seconds, which approximately corresponds to 1 km2. To derive the slope and 
aspect, the DEM data was processed within the R programming language (R Core Team, 2022). This comprehensive set of 
environmental parameters was chosen to capture the diverse range of factors that potentially influence the distribution of 
Rhododendron arboreum. 

Likewise, to address the issue of multicollinearity among the predictors, a Pearson correlation test was performed using the 
‘usdm’ package in R (Naimi, 2015). Variables with a correlation coefficient exceeding 0.7 and variation inflation factors (VIFs) 
greater than ten were identified and subsequently excluded from the analysis. Finally, 9 out of 23 variables were used to build the 
model: Mean Diurnal Range (BIO2), Temperature Seasonality (BIO4), Min Temperature of Coldest Month (BIO6), Annual 
Precipitation (BIO12), Precipitation of Driest Month (BIO14), Precipitation Seasonality (BIO15), Precipitation of Coldest Quarter 
(BIO19), Slope and Aspect (Table 1). 


Table 1 A complete list of predictor variables compiled for this study. Predictor variables used for the study are highlighted in bold 


text. 

Source Category Variable Unit 
BIO1 = Annual Mean Temperature °C 
BIO2 = Mean Diurnal Range (Mean of monthly (max temp-— 6 
min temp)) 
BIO3 = Isothermality (BIO2/BIO7) (100) % 
BIO4 = Temperature Seasonality (standard deviation x 100) i O 
BIO5 = Max Temperature of Warmest Month °C 
BIO6 = Min Temperature of Coldest Month 6 
BIO7 = Temperature Annual Range (BIO5-BIO6) °C 
BIO8 = Mean Temperature of Wettest Quarter =G 
BIO9 = Mean Temperature of Driest Quarter °C 

Worldclim | Bioclimatic BIO10 = Mean Temperature of Warmest Quarter °C 
BIO11 = Mean Temperature of Coldest Quarter °C 
BIO12 = Annual Precipitation mm 
BIO13 = Precipitation of Wettest Month mm 
BIO14 = Precipitation of Driest Month mm 
BIO15 = Precipitation Seasonality (Coefficient of Variation) % 
BIO16 = Precipitation of Wettest Quarter mm 
BIO17 = Precipitation of Driest Quarter mm 
BIO18 = Precipitation of Warmest Quarter mm 
BIO19 = Precipitation of Coldest Quarter mm 
Aspect Degree 

_ | Slope Degree 

SRTM Topographic - 

Elevation km 


Future Environmental Parameter 
To assess the response of Rhododendron arboreum to future climate conditions, we utilized the sixth version of the Model for 
Interdisciplinary Research on Climate (MIROC), also known as MIROC6. MIROC6 is a global circulation model (GCM) submitted 
for the Coupled Model Inter-comparison Project Phase 6 (CMIP6). MIROC's circulation has been demonstrated to effectively 
capture rainfall and temperature variability across South Asia, as highlighted by (Mishra et al., 2014). This characteristic makes 
MIROC particularly suitable for investigating the climate dynamics and patterns in the region under consideration. 

It is advisable to select multiple distinct future climate patterns to mitigate potential uncertainties and errors associated with 


future climate projections (Ncube et al., 2020). Hence, we chose bioclimatic variables representing two Green House Gas (GHG) 
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concentration paths based on Shared Socioeconomic Pathways (SSP): SSP2-4.5, which is based on an intermediate emission 
scenario, and SSP5-8.5, which is based on a very high emission scenario. By incorporating these different pathways, we aimed to 


account for a range of potential future climate scenarios and improve the reliability of our analysis. 


Species distribution modeling 

Species distribution modeling aims to predict spatial and temporal distribution of a species by analyzing the correlation between its 
occurrence or abundance and environmental conditions (Elith and Leathwick, 2009). Different tools, including statistical regression, 
machine learning, and geographic extrapolation, are used for this purpose (Shrestha and Shrestha, 2019). However, it is essential to 
note that the performance of these algorithms can vary significantly (Elith et al., 2010). Ensemble modeling, which involves 
simulating multiple sets of initial conditions, model classes, parameters, and boundary conditions, improves the reliability and 
accuracy of species distribution predictions (Araujo and New, 2007; Sagi and Rokach, 2018). Combining multiple approaches and 
accounting for algorithm uncertainties, ensemble modeling provides a consensus approach that outperforms individual models 
(Aguirre-Gutiérrez et al., 2017). In this research, we employed ensemble modeling to enhance the reliability and accuracy of species 
distribution predictions (Araujo and New, 2007; Thuiller et al., 2009). 

To conduct ensemble modeling, we utilized the 'biomod2' package within the R program (R Core Team, 2022). We selected ten 
algorithms: Multiple analysis regression splines (MARS), maximum entropy (MaxEnt), random forest (RF), boosted regression tree 
(BRT), generalized linear model (GLM), gradient boosting machine (GBM), classification tree analysis (CTA), artificial neural 
network (ANN), flexible discriminant analysis (FDA), surface range envelope (SRE) and generalized additive model (GAM) for 
creating the ensemble model. As actual absence data were unavailable, we randomly selected 10000 pseudo-absences for model 
training and evaluation. The models were calibrated using 80% of the occurrence points (presence and pseudo-absence) as training 
data, while the remaining 20% served as testing data to evaluate model performance. We performed seven evaluation repetitions 
using the ten selected model algorithms, resulting in 70 individual model runs. 

To evaluate predictive performance, we employed two widely accepted metrics: The area under the curve (AUC) of the receiver 
operating characteristic (ROC) plot and the true skill statistic (TSS) (Naimi et al., 2022). The AUC is a threshold-independent metric 
ranging from 0 to 1, categorizing models as excellent (AUC > 0.9), good (0.8 to 0.9), fair (0.7 to 0.8), or poor (0.6 to 0.7) (Ben-Rais- 
Lasram et al., 2010). The TSS, calculated as "sensitivity + specificity - 1", ranges from -1 to +1, with values below 0.4 indicating poor 
model discrimination (Allouche et al., 2006; Beaumont et al., 2016; Shrestha and Shrestha, 2019). Eight models (GLM, GBM, GAM, 
ANN, FDA, MARS, RF, MAXENT) with AUC values greater than 0.8 were selected for further analysis. 

The ensemble model was constructed using a committee averaging approach to combine predictions from those 8 individual 
models, aiming to enhance overall accuracy and reliability. We also determined the relative importance of each predictor variable 
used in the models. To determine the presence and absence based on the predicted probabilities, we employed the "Maximum 
Sensitivity plus Specificity" threshold, as described by (Liu et al., 2005; Naimi and Araujo, 2016). This threshold allowed us to 
convert the predicted probabilities into binary predictions. Subsequently, we calculated the areas of the expected current 
distribution of the species and the changes in these habitats between the current and future prediction scenarios. 


3, RESULTS 


Model performance 

The models showcased commendable performance based on their evaluation scores. The TSS value displayed some variation across 
the models, ranging from 0.63 for the FDA model to 0.80 for the GBM model. Similarly, the AUC value showed a range of scores, 
with the ANN model scoring 0.83 and the GBM model achieving an impressive score of 0.94. The final ensemble model, which 
incorporated eight carefully selected algorithms, outperformed the rest with an outstanding AUC value of 0.97 and a TSS value of 
0.77 (Table 2). These results genuinely highlight the strong performance of the ensemble model and its potential in predictive 
modeling. 


Table 2 Evaluation metrices of algorithms used in the study. 
Model | ANN | FDA | GAM | GBM | GLM | MARS | MAXENT | RF _ | Ensemble 
AUC | 0.83 0.88 | 0.91 0.94 | 0.89 | 0.90 0.93 1.00 | 0.97 
TSS 0.64 0.63 | 0.76 0.80 | 0.70 | 0.70 0.74 1.00 | 0.77 
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Variables contribution 

Two variables stood out as the most influential: the minimum temperature of the coldest month (BIO6) and the annual precipitation 
(BIO12), as shown in (Table 3). These two variables contributed significantly, accounting for 69.52% and 13.20% of the model's 
importance, respectively. Another important variable was the temperature seasonality (BIO4), which had a contribution of 9.55%. 
Interestingly, the bioclimatic variables played a more crucial role in determining the potential distribution of Rhododendron arboreum 


compared to the topographic variables. The cumulative contribution of the bioclimatic variables was remarkably high at 99.17%. 


Table 3 Importance and percentage contribution of predictor variables. 


Variable codes | Mean Variable Importance | Percentage 
BIO2 0.0178 1.67 

BIO4 0.1025 9.55 

BIO6 0.7465 69.52 
BIO12 0.1417 13.20 
BIO14 0.0069 0.64 

BIOI5 0.0186 1.73 

BIO19 0.0308 2.87 

Sope 0.0054 0.51 
Aspect 0.0034 0.32 


Current and Future Potential Distribution 

By analyzing various environmental variables and species occurrence data, we were able to generate a comprehensive 
understanding of the current distribution of Rhododendron within the country. The findings revealed substantial variations in 
suitability across Nepal, with a total of 115,317.1 km2 (78.17%) identified as unsuitable for the species, while 32,198.89 km2 (21.83%) 
were classified as suitable habitats, as shown in (Figure 2). The province-wise analysis of the current habitat distribution for 
Rhododendron arboreum revealed interesting patterns across Nepal. Province No. 1 exhibited a substantial suitable habitat area of 
6503.48 km2, accounting for 24.62% of its total area. This indicates a significant presence of favorable conditions for the species 
within this province. On the other hand, Province No. 2 did not exhibit any suitable habitat for Rhododendron arboreum, suggesting 
limited or unfavorable environmental conditions for the species in this region. Moving to Bagmati Province, we found a suitable 
habitat area of 6585.56 km2, representing approximately 32.11% of its total area. 

This province demonstrates a relatively high potential for supporting the growth and distribution of Rhododendron arboreum. 
Similarly, Gandaki Province showcased a suitable habitat area of 5064.26 km2, accounting for 22.88% of its total area. This suggests 
that the species can thrive in various parts of this province. Lumbini Province exhibited a suitable habitat area of 3663.27 km2, 
representing 18.89% of the province's total area. Although relatively more minor than other provinces, it still provides a significant 
suitable habitat for Rhododendron arboreum. Karnali Province showcased a suitable habitat area of 5044.67 km2, comprising 16.36% of 
its total area. Despite the lower percentage, the province harbors substantial areas where the species can flourish. Sudurpaschim 
Province demonstrated a notable suitable habitat area of 5337.66 km2, representing 26.45% of its total area. This indicates a 
relatively high potential for supporting Rhododendron arboreum within this province. 

The study examined the current distribution of Rhododendron in Nepal and made projections for its future suitable habitat 
under two different climate scenarios: SSP2-4.5 and SSP5-8.5. The analysis revealed that the species faces substantial changes in 
habitat suitability by the years 2050 and 2070, as shown in (Figure 3). 

Considering the current distribution patterns, which indicate that 115,317.1 km2 are unsuitable and 32,198.89 km2 are suitable 
for Rhododendron, we observed significant shifts in suitable and unsuitable areas. For the year 2050, under the SSP2-4.5 scenario, 
the models projected an increase in unsuitable areas, covering approximately 121,234.7 km2, and a decrease in suitable habitat, with 
26,281.3 km2 remaining. Conversely, under the SSP5-8.5 scenario for the same year, the unsuitable areas were projected to expand 
further, encompassing 123,858.3 km2, while the suitable habitat decreased to 23,657.7 km2. Looking ahead to 2070, under the SSP2- 
4.5 scenario, the models predicted a further decline in suitable habitat, with only 23,260.4 km2 remaining, and an increase in 
unsuitable areas, which expanded to 124,255.6 km2. 
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Figure 2 Predicted current distribution of Rhododendron arboreum in Nepal 
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Figure 3 Predicted future distribution of Rhododendron arboreum in Nepal. 


Species 24, e81s1590 (2023) 


7 of 14 


ANALYSIS ARTICLE OPEN ACCESS 


Similarly, under the SSP5-8.5 scenario, the projections indicated a significant reduction in suitable habitat to 19,913.73 km2, 
while the unsuitable areas expanded to cover 127,602.3 km2. These findings suggest a concerning trend of declining suitable 
habitats for Rhododendron in the future, with both climate scenarios pointing towards a potential decrease in suitable areas and an 
expansion of unsuitable regions. The projections indicate that the species may face increasing challenges to its survival and 
persistence, particularly under the SSP5-8.5 scenario. 


Habitat Gain and Loss 

We observed losses and gains in suitable habitats by comparing the current distribution patterns with the projected scenarios for 
the years 2050 and 2070 under the SSP2-4.5 and SSP5-8.5 climate scenarios (Figure 4). For the year 2050 under the SSP2-4.5 scenario, 
the models predicted a loss of approximately 10,710.39 km2 and a gain of 4,792.78 km2 in suitable habitats. Moreover, an area of 
21,488.5 km2 was projected to remain stable. In contrast, under the SSP5-8.5 scenario for the same year, the projections indicated the 
more extensive loss of 14,237.64 km2 and a smaller gain of 5,696.41 km2. Additionally, an area of 17,961.26 km2 was expected to 
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Figure 4 Habitat gain and loss of Rhododendron arboreum under different climate change scenarios for the future. 


Looking ahead to 2070, under the SSP2-4.5 scenario, the models predicted a further loss of 14,985.52 km2 and a gain of 6,047.03 
km2 in suitable habitat. Meanwhile, an area of 17,213.37 km2 was projected to remain stable. Similarly, under the SSP5-8.5 scenario, 
the projections indicated a more significant loss of 19,790.34 km2 and a more considerable gain of 7,505.18 km2. Additionally, an 
area of 12,408.56 km2 was expected to remain stable. 


Elevation range shift 

Our findings indicate that Rhododendron arboreum predominantly occupies hilly areas in Nepal within an elevation range of 1080 to 
3730 meters during the current period. Our projections for the year 2050, specifically under the SSP2-4.5 scenario, present a 
noteworthy upward shift in elevation for Rhododendron arboreum. The projected suitable habitat elevation range expands from 1373 
to 4029 meters, implying a potential movement of the species toward higher altitudes, likely in response to anticipated climate 
variations. The widening of the suitable habitat elevation range signifies the species' remarkable adaptability to evolving 
environmental conditions. 
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Similarly, under the SSP5-8.5 scenario for the year 2050, the suitable habitat elevation range extends even further, ranging from 
1525 to 4203 meters. This indicates a more pronounced upward shift compared to the SSP2-4.5 scenario, emphasizing the 
heightened impact of climate change on the species distribution patterns. These findings underscore the importance of considering 
various climate scenarios to comprehend the full spectrum of potential influences on Rhododendron arboreum. Looking ahead to 2070, 
our projections under the SSP2-4.5 scenario indicate a continuation of the upward trend in elevation shifts. The estimated suitable 
habitat elevation range extends from 1525 to 4208 meters, suggesting a continuous movement of the species toward higher 
altitudes. The projected expansion of the suitable habitat range further emphasizes the need for long-term conservation strategies 
that account for the anticipated changes in the species distribution patterns. 
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Figure 5 Elevation range of suitable habitat for Rhododendron arboreum for current and future period. In the x-axis, 1 represents the 
current range, 2 represents the range for 2050 ssp2-4.5, 3 represents the range for 2050 ssp5-8.5, 4 represents the range for 2070 ssp2- 
4.5 and 5 represents the range for 2070 ssp5-8.5 


Under the SSP5-8.5 scenario for the year 2070, the suitable habitat elevation range undergoes a significant expansion, ranging 
from 1864 to 4355 meters. This substantial upward shift, as shown in Figure 5, highlights the considerable displacement of 
Rhododendron arboreum towards higher elevations in response to changing climatic conditions. 


4. DISCUSSION 


To develop effective species recovery plans within their natural habitats, an important step involves modeling the present habitats 
and predicting the future range of species (Kumari et al., 2022). To ensure the conservation and sustainability of species habitats in 
the future, it is crucial to employ habitat distribution modeling, which enables us to assess dispersal range and the impact of climate 
change on future habitat suitability (Mushtaq et al., 2021; Zhong et al., 2021). This approach helps us understand the intricate 
relationship between species distribution range, extinction, and the critical determinants of geographical range area (Cardillo et al., 
2008; Purvis et al., 2000). In this context, the present study takes a pioneering approach by mapping the existing habitats of 
Rhododendron arboreum in Nepal and forecasting how climate change will influence the suitability of these habitats in the future 
using an ensemble approach. 

In our study, we employed the ‘biomod2’ package to model the potential distribution of Rhododendron arboreum in Nepal. The 
‘biomod2’ package offers a unique advantage by integrating multiple modeling algorithms and demonstrating high predictive 
accuracy. Notably, prior research by Wani et al., (2022) utilized ‘biomod2’ for predicting the habitat suitability of Dactylorhiza 
hatagiera in the Himalayan region under projected climate change. Their study achieved impressive evaluation metrics, with an 
AUC value of 0.93 and a TSS value of 0.82. Similarly, other studies by Farashi and Erfani, (2018), Kumari et al., (2022), Shrestha and 
Shrestha, (2019) have successfully employed ‘biomod2’ for ensemble distribution modeling of various flora and fauna, yielding 


notable predictive accuracy. In our investigation, we evaluated the performance of our models using commonly accepted 
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evaluation metrics. We obtained an AUC value of 0.97 and a TSS value of 0.77, considered acceptable and consistent with the 
aforementioned studies. These results highlight our modeling approach's robustness and reinforce our findings' reliability. 

The predicted current distribution of Rhododendron arboreum aligns remarkably well with the known occurrence records, 
providing a reliable baseline for understanding the species’ habitat preferences. This species favors hilly regions at moderate 
elevations Anand et al., (2021), which is consistent with its established ecological characteristics. In our study, we have meticulously 
compiled a comprehensive dataset comprising 160 occurrence records encompassing a wide range of elevations, from 1241 to 3281 
meters. Notably, an exceptionally high number of occurrences are observed within the elevation range of 1500 to 3000 meters, as 
depicted in (Figure 6). 


occurence 


1500 2000 2500 3000 
elevation 


Figure 6 Histogram of occurrence points with respect to elevation. 


Our models confirm that the species thrives within an elevation range of 1080 to 3730 meters, which corresponds closely with 
the findings of Veera et al., (2019) in the western Himalayas of India. Their research suggests that the distribution of Rhododendron 
arboreum typically occurs between 1500 and 3500 meters above sea level. Similarly, our study revealed that, based on the current 
climate conditions, the potentially suitable areas are primarily concentrated in the eastern and western parts, while the central 
region exhibits relatively fewer suitable areas. The province-wise analysis of the current habitat distribution for Rhododendron 
arboreum in Nepal provides valuable insights into the species' ecological preferences and distribution patterns. The findings reveal 
notable variations in suitable habitat availability across different provinces, emphasizing the importance of considering regional 
factors in conservation planning. 

Province No. 1 emerges as a significant hotspot for Rhododendron arboreum, covering nearly a quarter of the province's total area. 
This highlights the province's favorable environmental conditions, such as suitable elevation, climate, and vegetation, that support 
the growth and proliferation of the species. Identifying such a sizable habitat area in Province No. 1 underscores the need for 
targeted conservation efforts to preserve and protect this vital stronghold of Rhododendron arboreum. In contrast, Province No. 2 
lacks a suitable habitat for the species. The absence of favorable conditions in this region might be attributed to various factors, 
including unsuitable elevation ranges, land-use patterns, or climatic conditions, as described by (Sharma et al., 2020). 
Understanding the limitations or constraints that prevent the establishment of suitable habitats in Province No. 2 is crucial for 
informing conservation strategies in areas where Rhododendron arboreum currently faces challenges. Bagmati Province and Gandaki 
Province exhibit considerable suitable habitat areas, with 6585.56 km2 (32.11% of total area) and 5064.26 km2 (22.88% of total area), 
respectively. 

These findings indicate that these provinces provide favorable conditions for Rhododendron arboreum, potentially contributing to 
the species' ecological resilience and conservation. The presence of significant suitable habitat areas in these regions highlights their 
importance as key conservation areas and calls for targeted conservation initiatives to preserve the existing populations and 
promote their long-term viability. Lumbini Province, Karnali Province, and Sudurpaschim Province also demonstrate varying 
degrees of suitable habitat availability for Rhododendron arboreum, comprising 18.89%, 16.36%, and 26.45% of their total areas. 
Despite variations in the percentage of suitable habitat, these provinces play crucial roles in supporting the species distribution and 
maintaining its ecological balance within Nepal's diverse landscapes. 

Furthermore, when considering future climate scenarios, we found that the potentially suitable areas are expected to decrease 
overall, indicating a higher probability of habitat loss than habitat gain, as shown in (Table 4). Our findings indicated that across all 


projected future climate scenarios, there was a notable increase in suitable habitats in the higher latitude areas, especially in Karnali 
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and Sudurpaschim province, as highlighted in the study (Shifting Habitats, 2020). However, habitat loss was consistent throughout 
the entire country. These projected losses in the habitat of Rhododendron arboreum highlight the need for proactive conservation 


measures to ensure the preservation and protection of suitable habitats. 


Table 4 Area of suitable habitat of Rhododendron arboreum for current and future period as well as area gained and lost in future. 
Current 2050 2070 
ssp2-4.5 ssp5-8.5 ssp2-4.5 ssp5-8.5 


& Not suitable | 115317.11 | 121234.72 | 123858.33 | 124255.60 | 127602.27 
Z | Suitable 32198.89 | 26281.28 | 23657.67 | 23260.40 | 19913.73 
& Gain - 4792.78 5696.41 6047.03 7505.17 

o 

Z Loss - 10710.39 | 14237.64 | 14985.52 | 19790.34 


The two most influential variables were the minimum temperature of the coldest month (BIO6) and the annual precipitation 
(BIO12). The prominence of these two variables suggests that they play a crucial role in shaping suitable habitats for Rhododendron 
arboreum. The minimum temperature of the coldest month likely affects the plant's ability to tolerate low temperatures, while 
annual precipitation is a vital factor in determining the availability of water resources. Together, these variables appear to influence 
the potential distribution of Rhododendron arboreum strongly. Additionally, we found that the temperature seasonality (BIO4) had 
a notable contribution of 9.55%. This suggests that the variation in temperature throughout the year, rather than just extreme cold 
temperatures, also significantly determines the suitable habitats for Rhododendron arboreum. 

These findings highlight the importance of considering bioclimatic variables when studying the potential distribution of 
Rhododendron arboreum. By focusing on these factors, we can better understand the environmental conditions required for the plant's 
survival and identify suitable habitats for conservation and management efforts. It is worth noting that our study focused on a 
specific species in a particular region. Future research could explore the influence of these variables on other species or in different 
geographical contexts to gain a more comprehensive understanding of the factors influencing species distributions. Furthermore, 
incorporating additional variables, such as soil characteristics or land use patterns, may provide further insights into the ecological 
requirements and potential threats to Rhododendron arboreum. 


5. CONCLUSION 


In conclusion, this research provides valuable insights into the distribution patterns of Rhododendron arboreum in Nepal, highlighting 
its cultural and ecological significance. Through applying species distribution modeling techniques and utilizing an ensemble 
approach, we have gained a comprehensive understanding of the current distribution and projected future changes in suitable 
habitats for this iconic species. Our findings indicate that a significant portion of Nepal's land area is unsuitable for Rhododendron 
arboreum, emphasizing the need for conservation efforts to protect the remaining suitable habitat. The future projections under 
different climate scenarios reveal potential gains and losses in suitable habitats, pointing to the dynamic nature of the species 
distribution and the impact of climate change. The implications of this research call for proactive conservation measures to 
safeguard Rhododendron arboreum in the face of climate change. 

By understanding the potential challenges and predicting future habitat changes, we can develop informed conservation 
strategies to preserve the species and its cultural and ecological importance. The study contributes to the broader field of species 
distribution modeling by showcasing the effectiveness of ensemble approaches in generating reliable predictions. It is a foundation 
for further research and conservation efforts focused on the region’s Rhododendron arboreum and other vulnerable plant species. By 
protecting and managing suitable habitats, we can ensure the long-term survival of Rhododendron arboreum, preserving its beauty, 
cultural significance, and the intricate ecological balance it contributes to Nepal's diverse landscapes. 
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